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^the  wide-angle  time-domain  parabolic  equation  (TDPE),  which  is  the  inverse  Fourier 
transform  of  the  wide-angle  parabolic  equation  (PE),  is  derived.  A  numerical  solution  for  the 
model  is  described  and  a  benchmark  calculation  is  presented.  The  narrow-angle  TDPE  is  also 
considered  and  its  error  is  analyzed  and  compared  with  the  error  of  the  narrow-angle  PE.  The 
TDPE  is  compared  with  the  progressive  wave  equation,  which  is  shown  to  be  restricted  to 
narrow-angle  propagation  for  practical  purposes.  In  the  sediment,  attenuation  is  assumed  to 
depend  linearly  on  frequency  and  the  corresponding  causal  dispersion  law  is  assumed.  The 
model  is  used  to  show  that  the  effect  of  sediment  dispersion  on  pulse  propagation  in  the  ocean 


can  be  significant.  '*■  -  -  Z-  ■ 

PACS  numbers:  43.30.Bp,  43.30.Ma,  43.20.Bi 


INTRODUCTION 

Time-domain  approaches  are  useful  for  modeling 
broadband  acoustic  propagation.  For  example,  suppose  that 
a  sequence  of  snapshots  of  the  acoustic  pressure  in  two  spa¬ 
tial  variables  is  desired  to  study  the  evolution  of  a  pulse  in 
time.  With  frequency-domain  approaches  to  this  problem,  it 
is  necessary  to:  (a)  Fourier  decompose  the  source  function; 
(b)  solve  the  propagation  problem  for  each  frequency;  (c) 
store  and  manage  the  solution  for  each  frequency;  and  (d) 
perform  a  three-way  sum  over  frequency  and  space  for  each 
snapshot.  Errors  due  to  approximations  and  round-off  occur 
in  steps  ( a ) ,  ( b ) ,  and  ( d ) .  Step  ( a )  requires  the  selection  of  a 
frequency  spacing  and  a  frequency  band.  Step  (b)  requires 
the  selection  of  grid  spacings  for  each  frequency.  Step  (c) 
can  be  difficult  for  broadband  problems.  Even  if  step  (b) 
requires  less  computer  time  (CPU)  than  a  time-domain  cal¬ 
culation,  it  is  possible  that  step  (d)  will  offset  the  advantage. 
For  an  analogous  problem, 1  inverting  a  modal  decomposi¬ 
tion  required  several  times  as  much  CPU  as  performing  the 
calculation  for  each  of  the  modes. 

In  shallow  water,  it  is  essential  for  an  underwater  pulse 
propagation  model  to  handle  bottom  interaction,  range-de¬ 
pendence,  and  wide-angle  propagation.  Since  the  ocean  is  an 
inhomogeneous  waveguide  of  variable  depth,  a  realistic 
propagation  model  must  handle  both  bottom  interaction 
and  range  dependence.  Benchmark  studies2  indicate  that 
wide-angle  capability  is  important  in  underwater  acoustic 
modeling.  The  parabolic  equation1  (PE)  method  is  ideally 
suited  to  handle  range  dependence.  With  the  development  of 
wide-angle  capability, the  PE  method  became  the  most 
useful  frequency-domain  tool  available  for  bottom-interact¬ 
ing  propagation.  Two  time-domain  methods  related  to  the 
PE  method  have  been  developed.  The  progressive  wave 
equation^*  (PWE)  advances  the  acoustic  pressure  in  time. 
The  PWE  has  been  extended  to  handle  nonlinear  propaga¬ 
tion7  and  bottom  interaction.1'  The  time-domain  parabolic 
equationh'y  (TDPE)  is  the  inverse  Fourier  transform  of  the 
PE.  It  advances  the  acoustic  pressure  in  range.  More  effort 


;  ,  I  if  ....  /,>  *r.-  ■ 

has  gone  into  the  development  of  the  PWE,  perhaps  because 
it  is  more  natural  to  march  a  solution  of  the  wave  equation  in 
time  rather  than  range. 

In  this  article,  the  TDPE  is  extended  to  handle  wide- 
angle  propagation,  and  an  energy  argument  is  presented  that 
shows  that  the  PWE  is  not  useful  for  wide-angle  propaga¬ 
tion.  A  benchmark  calculation  is  presented  and  the  error  of 
the  narrow-angle  TDPE  is  compared  with  the  error  of  the 
narrow-angle  PE.  Sediment  attenuation  is  assumed  to  de¬ 
pend  linearly  on  frequency,  which  agrees  with  experimental 
results  involving  various  materials  and  frequencies. 10-12  The 
corresponding  causal  sediment  dispersion  relation,"  which 
has  been  validated  experimentally,"  is  also  assumed  in  the 
sediment.  A  calculation  is  presented  to  demonstrate  that 
sediment  dispersion  can  have  a  significant  effect  on  pulse 
propagation  in  the  ocean. 


I.  THE  FREQUENCY-DOMAIN  PARABOLIC  EQUATION 

A  time-harmonic  steady  state  is  assumed  and  the  acous¬ 
tic  pressure  p  is  factored  as  p(\,t )  =/,(x)exp(  —icot), 
where  t  is  time,  x  is  the  Cartesian  position  vector,  and  to  is  the 
circular  frequency.  The  complex  pressure  P  is  assumed  to 
satisfy  the  pressure-release  boundary  condition  P  =  0  at  the 
ocean  surface,  the  outgoing  radiation  condition  at  infinity, 
and  the  reduced  wave  equation14 

pV'[(l/p)VP  ]  +  K  :P  =  -  4trd(x  -  x,t).  (1) 

where  the  point  x„  is  t  he  source  location.  The  complex  wave- 
number  K  =  k  +  icrP  |k  ]  is  used  to  account  for  sediment  loss. 
Absolute  value  is  used  so  that  energy  loss  occurs  in  the  direc¬ 
tion  of  propagation.  The  wavenun  ber  is  k  =  to/c, 
o=  (40/7  log,,,  e)  ',  /S  is  the  attenuati  >n  in  decibels  per 
wavelength  (dB/i),  p  is  the  density,  and  c  is  the  sound 
speed.  The  variable  density  term  is  due  to  Bergmann. " 

To  reduce  to  two  spatial  dimensions,  we  assume  that 
azimuthal  variations  are  negligible.  Since  the  ocean  is  a 
waveguide,  energy  propagating  from  a  source  exhibits  cylin- 
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drical  spreading.  Thus  it  is  often  beneficial  to  sol  ve  Eq.  ( 1 )  in 
cylindrical  coordinates  with  z  being  the  depth  below  the 
ocean  surface  and  r  being  the  horizontal  distance  from  a 
source  at  the  depth  z0.  Variations  in  range  are  assumed  to  be 
sufficiently  weak  so  that  dp/ dr  can  be  ignored,  which  simpli¬ 
fies  Eq.  ( 1 )  to 

cT-P  1  dp  dP  |  d2P  |  1  <?/>  !  K2p 

dr  p  dz  dz  dr 2  r  dr 

=  -  (2/r)  <5(r)<5(z-z0).  (2) 

The  propagation  angle  of  a  ray  is  defined  to  be  the  angle 
it  makes  with  the  ocean  surface.  We  define  aM  to  be  the 
maximum  angle  of  propagation  in  the  ocean.16  Rays  that 
propagate  with  angles  greater  than  aM  are  not  trapped  in  the 
ocean  and  thus  do  not  contribute  to  the  farfieid.  Since  the 
discontinuity  in  sound  speed  at  the  ocean  bottom  is  small, 
e  =  tan2  aM  <  1  and  k0  '|fc  -  £0|  =  0(e),  where  k„  is  an 
average  wavenumber  in  the  water  column.  We  define 
Q  =  r'l2P,  and  Eq.  (2)  becomes 


?!Q_±4LW  +  d^Q+Q_  +  K,Q  =  o. 

dr  p  dz  dz  dr 2  \r 1 


(3) 


We  assume  that  r  >  r0,  where  k0r0  >  1 ,  and  drop  the  0(r  2 ) 
term  in  Eq.  ( 3 )  to  obtain 

_ 

P 


dr 


±*PL&  +  *!e  +  K2Q  =  0. 

dz  dz  dr 2 


(4) 


We  solve  Eq.  (4)  for  d 2 /dr  2  in  operator  notation  and  take 
the  square  root  to  obtain  the  outgoing  operator 


d_ 

dr 

where 


/  ^ 

=  <*°  /!+  — 


+  i- 


L  = 


il 

dr 


\__dp_d_ 
p  dz  dz 


(5) 


(6) 


We  consider  a  plane  wave  traveling  with  the  vertical  angle 
n<aM: 

Q„  =  exp[/'£0(rcos  a  +  z  sin  nr)  ].  (7) 

Since dQ„ /dz  =  0(el/2)  andp depends  weakly  on  depth,  we 
are  motivated  to  assume  that 


Q=Q'(ze'n),  (8) 

p  =  p'(ze‘/:),  (9) 

where  Q'  and  p  are  independent  of  e.  Thus  k  „  2L  —  0(e) 
and  we  may  replace  the  square  root  in  Eq.  { 5 )  with  its  Taylor 
series  approximation  to  obtain  the  narrow-angle  PE  opera¬ 
tor 


d_ 

dr 


=  ikn  ^ 


1  + 


K2  -  k2,  +  L 
2  k2 


) 


(10) 


Replacing  the  square  root  with  its  Pade  approximation,  we 
obtain  the  wide-angle  PE  operator 

2(K 2  —  i 

IT; 


dr  V  3  k}t+K2  +  L) 


(11) 


We  define  O' by 


Q(r,z)  =  U(r,z)exp(ik0r).  (12) 

To  simplify  Eqs.  (10)  and  ( 1 1 ) ,  we  assume  that  cr/34  e.  For 
the  wide-angle  operator,  we  also  require  that 

k0  '| k  —  A()|  <ein  the  water.  This  assumption,  which  is  val¬ 
id  in  shallow  water,  leads  to  a  TDPE  that  is  easy  to  solve 
numerically.  Substituting  Eq.  ( 12)  into  Eqs.  ( 10)  and  (11), 
we  obtain  the  narrow-angle  PE 

i£  =  i(k-k0)U-aP\k()\U+~U  (13) 

dr  2  kt) 

and  the  wide-angle  PE 

iH  =  i(k-k0)0--oP\kl)\U+-^^-U.  (14) 

dr  4  k  2  +  L 

Since  an  outgoing  signal  is  trapped  in  the  water  column,  Eq. 
( 14)  should  be  valid  for  wide-angle  propagation  in  shallow 
water. 

We  solve  Eq.  (14)  with  the  method  of  alternating  direc¬ 
tions,17  which  requires  numerical  solutions  for  each  of  the 
following: 


dU_ 

dr 


-  i(k  -  k„)0. 


(15) 


dU 

dr 


-  ol3\k„\U, 


(16) 


j  dU_  d"U  _  J_  dp^  d  U 
°  dr  dr  dr  p  dz  dr  dz 


=  2  ik 


d2U 
o  ~rr 
dr 


-  2  ik 


±dp_dU 
p  dz  dz 


(17) 


Equations  (15)  and  ( 16)  can  be  solved  exactly.  We  apply 
Galerkin’s  method  to  reduce  Eq.  (17)  to 


R  — —  +  SU  =  0,  (18) 

dr 

where  R  and  S'  are  tridiagonal  matrices  and  U  is  the  vector 
containing  the  values  of  U  at  the  depth  grid  points.  This 
approach  is  effective  for  handling  piecewise  continuous  var¬ 
iations  in  p  and  k.'K  Details  regarding  the  entries  of  the  ma¬ 
trices  are  given  in  the  Appendix  and  in  Ref.  8.  Crank-Nicol- 
son  integration  is  used  to  solve  Eq.  (18);  Eq.  (13)  is  solved 
in  a  similar  fashion. 

We  demonstrate  the  validity  of  Eq.  ( 14)  and  the  nu¬ 
merical  method  with  a  benchmark  problem 1,1  for  which  data 
appear  in  Table  I.  Subscripts  w  and  b  stand  for  water  and 
bottom  and  zr  is  the  receiver  depth.  The  ocean  depth 
d  —  100  m  is  constant.  The  deep  source  excites  wide-angle 
modes.  Transmission  loss  data  obtained  with  a  normal 
modes  calculation  and  with  Eqs.  (13)  and  (14)  using  a 
Gaussian  PE  starter'  appear  in  Fig.  1.  The  results  of  Eq. 


TABLE  I.  Data  for  the  wide-angle  PE  benchmark  problem  CPU  for  each 
PE  calculation 


z„  —  99.5  m 

z,  ---  9<).5  m 

c„  =  1500  m/s 

rK  —  1500m/s 

p„  --  1  g/cm' 

n„  ~  o 

c„  =  1 590  m/s 

fi,.  =  1.2  g/cm’ 

/?,.  =  0.5  dB//i 

d  -  100  m 

f()  5(X)  <r  s  1 

Ar  -  1  m 

Ar  =  0.2  m 

ew  =  200  m 

CPU  -  4  5  mm 
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FIG.  1.  Wide-angle  PF.  benchmark.  Solid  curve  is  (a)  narrow -angle  PE 
result  and  (b)  wide-angle  PE  result.  Dashed  curve  is  normal  modes  result 


(13)  exhibit  phase  delay  errors  due  to  wide-angle  propaga¬ 
tion  The  results  of  Eq.  (14)  agree  well  with  the  normal 
modes  result,  which  shows  that  Eq.  (14)  is  a  valid  wide- 
angle  PE.  All  calculations  were  done  on  a  Digital  VAX-8650 
computer. 

II.  THE  TIME-DOMAIN  PARABOLIC  EQUATION 


inverse  transform  of  the  operator  —  crfi  ,w\/c„by  verifying  it 
for  a  single  frequency.  We  rewrite  Eq.  (14)  as 


*■) 

or  \  c  c„/ 


o/3\co\ 


2ia>c„  L 


U  _  '  '  i  u  _j_  fj 

Co  4ru"  -f-  C'0L 

(23) 

and  invert  the  Fourier  transform  to  obtain  the  wide-angle 
TDPE 


du 


Vc„  c  /  dt  1 TC(>  J 


u(t')  —  u(t) 

(r'~n2  ' 


dt' 


+  ■ 


lct)(d/dt)L 


(24) 


4  (d'Zdt-)-clL 

As  for  the  nonlinear  PWE  of  Ref.  7,  each  of  the  terms  on  the 
right  side  of  Eq.  (24)  accounts  for  a  specific  physical  pro¬ 
cess;  and  they  are  referred  to  ( from  left  to  right )  as  the  re¬ 
fraction  term,  the  attenuation  term,  and  the  wide-angle  dif¬ 
fraction/density  term. 

The  TDPE  can  be  initialized  at  r  =  r„  by  the  homoge¬ 
neous  half-space  field"' 

p.(w), 2-/(, -^)  -i-).  (25, 

d  \  =  rl  +  (z±z„):,  (26) 

where/ (t)  is  the  source  function.  The  half-space  field  satis¬ 
fies  the  pressure  release  boundary  condition  at  the  ocean 
surface,  and  it  accounts  for  the  direct  arrival  and  the  surface- 
reflected  arrival  without  accounting  for  refraction,  loss,  or 
bottom  reflections.  This  starter  is  accurate  because  refrac¬ 
tion  and  attenuation  are  weak  and  can  be  neglected  near  the 
source.  Furthermore,  rays  that  reflect  from  the  ocean  bot¬ 
tom  near  the  source  propagate  at  large  angles  Thus  they  are 
not  trapped  in  the  oceanic  waveguide  and  do  not  affect  the 
farfield. 

We  define  q(r,z,t)  by  q  =  r'  "2p  and  observe  that 


q(r,z.t ) 


-r. 


Q(r,z,b))e*  p(  —  iiot)d(o, 


Q(r,z,u>)  - 


f. 


q{r,z,t)exp(icot)dt 


(27) 


(28) 


We  define  u  by 

X 

u(r,z,t)  =  j  U(r,z,(0)ex  p(  —  i(ot)dw. 

(19) 

u{r,z.t)  =  Q{  r,z,  ot)exp 

-  m  (,  +  -) 

U(r,i 


1  r 

\z,o) )  = -  (  u  ( r,z, 

2  y  J  „ 


t)exp(icoi)dt. 


We  define  c„  =  (o/km  rewrite  Eq.  (13)  as 
dU 
dr 


-v(i-i) 

V  c  cj 


U 


M  L,  '// 

c,  2(o 


(20) 


V,  (21) 


and  invert  the  Fourier  transform  to  obtain  the  narrow-angle 

TDPE 


du  f  1 _ 1  A  du 

dr  Vc„  c  /  dt 

■r 


+ 


£0 

irc„ 

u( t 1 )  -  u(/) 


,  ,  c„  L 
dt  +  —z—  u. 


(22) 


(/'-/)-  2  d/dt 

The  integral  in  Eq.  ( 22 )  exists  as  the  Cauchy  principal  value. 
One  can  show  that  the  integral  operator  in  Eq.  (22)  is  the 


dco. 

(29) 

From  Eq.  (29),  we  deduce  that  u{r,z,t)  ~  q(r,z,t  +  r/c,,). 
Thus  u  has  the  values  of  q  in  a  reference  frame  that  moves  in 
time  and  tracks  an  outgoing  signal.  In  contrast,  the  PWE 
involves  v(r,z,t)  =  q(r  +  cltt,z,t),  where  r  has  the  values  of  q 
in  a  reference  frame  that  moves  in  space. 

To  solve  the  TDPE  numerically,  the  source  function 
/(f)  is  assumed  to  have  compact  support.  A  time  window 
f,  <  t  <  f2  that  contains  the  signal  at  all  times  is  chosen.  This 
is  possible  since  the  outgoing  signal  is  tracked  by  the  time 
window  The  boundary  condition  u  —  0  is  used  at  the  pres¬ 
sure  release  surface,  deep  within  the  sediment  at  z  —  zM 
( from  which  no  energy  returns  to  the  water  column  due  to 
attenuation),  and  after  the  signal  has  passed  the  observer  at 
i  —  f:.  The  boundary  conditions  u  —  du/dr  ~  0  are  used  be¬ 
fore  the  signal  is  detected  at  f  = 
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We  solve  the  wide-angle  TDPE  with  the  alternating  di¬ 
rections  method  with  the  splitting  used  to  solve  the  wide- 
angle  PE: 


du 

+  (± 

i  ' 

)^  = 

0, 

(30) 

dr 

\  c 

ri>. 

/  dt 

du 

_  Ojg 

r 

u(t’)  - 

u(t) 

dt  ’  y 

(31) 

dr 

irc0 

j  * 

(f- 

t y- 

4 

d}u 

d 

*u  1 

dp 

d2u 

■» 

Co 

dr  dt 1 

dr  dr  p 

dz 

dr  dz 

_  2 

d'u 

2 

dp 

d2u 

(32) 

dz2  dt 

CoP 

dz 

dz  dt 

The  Lax-Wendroff  method10  is  used  to  solve  Eq.  (30), 
which  is  a  first-order  hyperbolic  equation. 

Since  Eq.  ( 16)  can  be  solved  exactly,  we  define  the  range 
increment  Ar  and  p  =  <7/3  Ar/c„  and  solve  Eq.  (31 )  as  fol¬ 
lows: 


U(r  +  Ar.z,<o)  =  U(r,z,(o )exp(  —  ym)i ). 
u(r- 1-  A r,z,t) 


U{r,z.  «)exp(  —  yj<u)  )exp(  —  ia>t)dio. 


(33) 


(34) 


Substituting  Eq.  (20)  into  Eq.  (34)  and  interchanging  the 
order  of  integration,  we  obtain 

u(r+Ar,z./)=i  f  -  u(r'z:' ’>  dt  \  (35) 

it  J,,  y-  +  (t  ‘  -  t)- 

The  kernel  of  the  integral  operator  of  Eq.  (35)  is  large  near 
t '  =  t  because  it  converges  to  <5  ( / '  —  t )  as  Ar — 0.  This  singu¬ 
lar  behavior  makes  the  operator  local  and  thus  numerically 
efficient.  Since  the  kernel  goes  to  zero  rapidly  away  from 
t'  —  t,  the  integration  limits  can  be  collapsed  to  a  small  inter¬ 
val  containing  t '  =  I.  We  approximate  the  integral  over 
(/,„— A//2,  r,„+A//2)  by  replacing  u(r,z,t)  with 
u(r,z,tm  )  and  integrating  the  kernel  exactly  to  obtain 


u(r+  \r,z,t„ ) 


tan 


—  tan  1 


u(r,z,tm ), 

(36) 


where  the  time  grid  points  are  t„  =  n  Ar  and  A t  is  the  time 
increment. 

Galerkin’s  method  is  used  to  reduce  Eq.  (32)  to 


Air  +  B~fh 

dr  dr  dt  ‘ 


+  C— =  0, 

dt 


(37) 


where  A ,  B,  and  C are  tridiagonal  matrices  and  u  is  the  vector 
containing  the  values  of  u  at  the  depth  grid  points.  Details 
regarding  the  entries  of  these  matrices  are  discussed  in  the 
Appendix.  Crank-Nicolson  integration  is  used  to  advance 
Eq.  (37)  in  range  using  centered  differences  in  t.  Since  the 
energy  flow  due  to  geometric  dispersion  is  from  t,  to  t2,  it  is 
necessary  to  sweep  from  f ,  to  t2. 

We  demonstrate  the  validity  of  Eqs.  (22)  and  (24)  with 
a  benchmark."  We  use  the  data  from  Table  II  with  the  Gaus¬ 
sian  source  function 


TABLE  II.  Data  for  the  wide-angle  TDPE  benchmark  problem  CPU  for 
wide-angle  TDPF  calculation. 


z0  =  75  m 
cu.  =  1 500  m/s 
cb  =  1600  m/s 
d  =  d{r) 

Ar  =  5m 
to  =  IOOtt  s  ' 


z,  =  25  m 
=  1  g/cm’ 
p,.  =  1.5  g/cnv 
150  s  1 
Az  =  2  m 
zM  =  300  m 


c„  =  1 500  m/s 

P„=  0 
0„  =  0.5  dB/Z 
ru  —  50  m 
At  =  2/3  ms 
CPU  =  3.75  h 


/(t)  =  exp[  -  (v/)2],  (38) 

where  v  =  1 50  s _  1 .  The  ocean  depth  is  200  m  for  r  <  4  km. 
linearly  sloping  from  200  to  50  m  over  4  km  <  r  <  8  km,  and 
50  m  for  r>  8  km.  To  prevent  reflections,  a  layer  of  sediment 
lOOAz  thick  is  added  below  z  =  zM  in  which  P  increases 
linearly  to  10  dB/A.  We  have  found  this  approach  effective 
and  use  it  in  all  calculations.  The  results  of  Eq.  (24)  appear 
in  Fig.  2.  To  obtain  a  sequence  of  snapshots  of  the  acoustic 
pressure,  we  convert  the  horizontal  axis  from  time  to  range 
using  the  approximation 

p(r  +  8r,z,t)  == p(r,z,t  —  <5r/c0),  (39) 

which  is  valid  for  small  propagation  angles  and  small  dr. 
With  this  conversion,  it  is  easier  to  describe  the  snapshots. 
The  solid  contours  correspond  top  >  0;  the  dashed  contours 
correspond  to  p  <  0. 

In  the  water  column,  the  field  consists  of  a  sequence  of 
fronts  involving  multiple  reflections  from  the  ocean  surface 
and  bottom.  Since  the  reflection  coefficient  at  the  ocean  sur¬ 
face  is  -  1  and  the  reflection  coefficient  at  the  ocean  bottom 
is  approximately  —  1  for  small-angle  incidence,  the  multi¬ 
ply  reflected  arrivals  alternate  between  solid  and  dashed 
contours.  Energy  flows  to  the  left  due  to  geometric  disper¬ 
sion.  As  time  increases,  the  fronts  are  squeezed  together. 
Large  amounts  of  energy  penetrate  into  the  ocean  bottom  in 
the  upsfope  region. 

The  error  of  the  narrow-angle  TDPE  is  analogous  to  the 
error  of  the  narrow-angle  PE  as  both  involve  delays.  Distinct 
features  in  the  narrow-angle  transmission  loss  curve  in  Fig.  1 
appear  at  larger  ranges  than  they  should.  For  example,  the 
large  null  that  occurs  before  r  —  7  km  appears  well  beyond 
r  =  7  km.  It  is  evident  from  the  waveforms  appearing  in  Fig. 
3  that  energy  is  dispersed  too  fast  and  squeezing  is  delayed 
for  the  narrow-angle  TDPE.  The  errors  are  small  for  the  first 
arrivals.  However,  the  agreement  gets  progressively  worse 
for  the  later  arrivals  since  the  propagation  angle  increases 
with  arrival  time.  The  agreement  improves  as  wide-angle 
energy  is  cut  off  with  decreasing  depth. 

To  obtain  a  TDPE  benchmark,  we  approximate  a  time- 
harmonic  source  by 

sin(tof)  a"f  (/  —  t "),  (40) 

where  <u  =  100 tt  s  V  By  superposition,  the  time-harmonic 
response  pcw  is  approximated  by 

p^(r,z,t)»'£a"p(r,z,t  -  f),  (41) 

where p  is  the  response  to  /  Details  regarding  the  constants 
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a "  and  t"  and  the  approximation  of  a  function  by  Gaussians 
are  given  in  Ref.  8.  Transmission  loss  data  obtained  with  Eq. 
( 14)  as  well  as  Eqs.  (22)  and  (24)  using  Eq.  (41 )  appear  in 
Fig.  4.  Phase  errors  are  evident  for  the  narrow-angle  TDPE 
calculation.  As  for  the  waveforms  in  Fig.  3,  the  errors  de¬ 
crease  as  wide-angle  propagation  is  cut  off.  The  agreement  is 
good  for  the  wide-angle  TDPE  calculation. 

III.  COMPARISON  OF  THE  TDPE  AND  THE  PWE 

The  TDPE  and  the  PWE  are  special  cases  of  a  one- 
parameter  family  of  methods  for  solving  the  wave  equation. 
For  c  =  p  =  1  and  (3  =  0,  the  inverse  Fourier  transform  of 
Eq.  (4)  is 

+  iis.  ,=  .  ,42) 

Oz  3r  ■  3l 


7750  7850  7950  8050 


RANGE(m) 


9750  9850  9950  10050 


RANGE(m) 


FIG.  2.  Snapshots  of  the  pressure  field  obtained  with  wide-angle  TDPE 
near  (a)  r  =  2  km.  (b)  r  —  4  km,  (c)  r  =  6  km,  (d)  r  =  8  km,  and  (e) 
r=  10  km.  Horizontal  line  is  the  ocean  bottom  (gradual  slope  is  not 
shown ). 


In  the  derivation  of  the  TDPE.  we  defined  a  new  dependent 
variable.  In  the  following  analysis,  it  is  convenient  to  define 
new  independent  variables.  We  introduce  <  =  (t  —  r)/ 2, 
//  —  r  cos  <b  -f 1  sin  <b.  and  C  =  z,  which  transform  Eq.  (42 ) 
into 

d  q  ,  ,  .  ,  3  q 

- -  —  ( cos  ©  4-  sin  ©)  - — 

3C '  <hj 

-f  (cos’  (b  -  sin  <t>)  Q .  —  0.  (43) 

<hr 

We  require  that  q  vanish  at  c  —  c„,  C  -  and  C  =  C,  and 
that  3q/<X'  vanish  at  c  -  2',,.  At  >/  j/,„  we  impose  the  initial 

condition  q  q„.  I  he  geometry  is  illustrated  in  Fig.  5. 

Since  the  outgoing  field  is  tracked,  and  it  propagates 

Si 
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FIG.  A  Waveform  detected  by  received  at  ;  =  25  m  and  (a)  r  -  2  km.  (b) 
r  -  4  km.  (c)  r  -  6  km.  (d)  r  =  8  km,  and  (e)  r  •-  10  km.  Solid  curve  is 
wide-angle  TLl’F.  result.  Dashed  curve  is  narrow-angle  TDPE  result 


with  small  vertical  angles,  the  dominant  operator  in  Eq. 
( 43 )  is  9  /<?£.  Thus  we  assume  that  9  /dr/  4,9  /J£  4 9  / 9 and 
drop  the  third  term  to  obtain 

3  2  32 

— %  —  (cos  <t>  +  sin  <b) - —  =  0.  (44) 

9C,  ‘  9£  9t] 

The  narrow-angle  TDPE  corresponds  to  <t>  =  0  in  Eq.  (44). 
The  narrow-angle  PWE  corresponds  to  <t>  =  it/ 2.  Thus  the 
narrow-angle  TDPE  and  the  narrow-angle  PWE  have  the 
same  canonical  form. 

Differentiating  Eq.  (43)  with  respect  to  £  and  Eq.  (44) 
with  respect  to  r],  we  obtain 


9  'q  .  ,  .  d'q 

- 2—  —  (cos  d>  +  sm  0)  — 

9£  9c,'  9c'  9i] 


+  (cos:  <t>  —  sin:  <b) 


9'g 


9c  9r) 


T  “  0, 


(45) 


3  (  3  1 

- - - (cos  <f>  4-  sin  Ml  —  ?  ^  0. 

Jr/  9c  1  9§  9  V 


(46) 


Wesoiveforr?  'q/9c  9rf  in  Eq.  (46)  and  substitute  the  result 
into  Eq.  (45)  to  obtain 
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FIG  4  Wide-angle  TDFE  benchmark.  Solid  curve  is  (a I  narrow-angle 
TDPE  result  and  (b)  wide-angle  TDPE  result.  Dashed  curve  is  wide-angle 
PE  result. 


d'g  .  ,  .  ..  d'g 

- -  -  (cos  d  -f  sin  0) - — 

d£  dC  -  d£ 2  Dr, 

D  .3 

+  (cos  <b  —  sin  <t>) - -  =  0.  (47) 

drid£- 

As  the  sign  of  the  third  term  changes  at  d  =  tt/ 4,  the  canoni¬ 
cal  form  of  Eq.  ( 47 )  changes. 

The  wide-angle  TDPE  is  obtained  by  taking  <f>  =  0  in 
Eq.  (47) 


FIG  5.  Geometry  of  coordinate  systems.  Signal  propagates  within  the  diag¬ 
onal  strip. 


d  *q  d  'q 


+  -I'g-  -=  o. 

dc,  ~  dr)  dr)  d£  2 


(48) 


d£dC2 

We  investigate  the  stability  of  Eq.  (48)  with  the  energy 
method.’1  Due  to  the  boundary  conditions  at  c  =  t(>,  Eq. 
(48)  is  equivalent  to 

A=  f  fJX;#-#'.  (49) 

Dr,  Je„  d£  2  ~  J4l.  Jf„  dr,  d;  ’  -  - 

We  multiply  Eq.  (49)  by  q  and  integrate  over  f  using  inte¬ 
gration  by  parts  to  obtain 


££(#)’ 


a 

dr) 


=  -2 
We  define  the  energy  E  by 


d§  "  d£ '  +  q2 


dC 


OL,  . 


E^,-CHSAf)dr,ir+'1 


(50) 


di.  (51) 


Since  DE  /dr)  <  0  hv  Eq  ( 50) ,  we  deduce  that  the  wide-angle 
TDPE  is  well-posed.  A  similar  argument  holds  for  6  <  7t/4. 
Thus  Eq.  (47)  is  well-posed  for  4>  <  ;r/4. 

The  wide-angle  PWE  is  obtained  by  taking  <b  =  n/2  in 
Eq.  (47) 

d  2q  d  'q  d  <q 


d£,  dt,  2  df-dr)  dr)d$ 2  °’ 


dr, 


d£ 


_iu  (#)'**'• 

We  define  the  energy  E  =  |F|,  where 


F(Z,V)  = 


:£(#)■ 


dr  df-q2 


(52) 


(53) 


d$.  (54) 


If  F>  0  initially,  dE/dr)>0  by  Eq.  (53).  Thus  the  wide- 
angle  PWE  is  ill-posed  as  is  Eq.  (47)  for  6  >  rr/4. 

To  understand  why  the  wide-angle  PWE  is  ill-posed,  we 
replace  the  initial  data  at  g  =  with  the  boundary  data 
q  =  Oat  f  =  £n  and  We  multiply  Eq.  (52)  by  q  and  inte¬ 
grate  over  both  C  and  £  using  integration  by  parts  to  obtain 

<»• 

Thus  the  wide-angle  PWE  is  conservative  and  well-posed  as 
a  two-point  boundary  value  problem  for  <p  >  7r/4.  With  these 
boundary  conditions,  the  PWE  counterpart  of  Eq.  ( 32 )  can¬ 
not  be  solved  by  sweeping  in  t  making  it  impractical  for  nu¬ 
merical  calculations. 

We  now  illustrate  a  problem  associated  with  the  fact 
that  the  PWE  is  an  initial  value  problem  in  time  as  opposed 
to  range.  The  maximum  range  at  which  ph  may  be  applied 
decreases  as  the  source  gets  closer  to  the  bottom.1'’  Thus  ph 
cannot  be  used  to  initialize  the  field  over  a  wide-range  win¬ 
dow  for  a  deep  source.  We  modify  the  problem  described  in 
Table  II  by  moving  the  source  down  to  z  =  180  m  and  using 
the  half-space  field 


Ph(r,z)  =  (1  /d  )exp (ikud  ) 


( \/d  j  )exp (ik0d ,  ) 
(56) 
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as  an  initial  condition  at  r  =  r0.  Transmission  loss  obtained 
with  Eq.  (14)  appears  in  Fig.  6forr„  =  50  and  250  m  and  for 
a  Gaussian  PE  starter.  The  error  for  the  case  ru  =  250  m 
shows  that  ?h  is  not  valid  to  250  m.  Based  on  the  waveforms 
in  Fig.  2,  which  are  approximately  250  m  in  width,  a  window 
of  this  width  is  reasonable.  Thus  it  would  be  difficult  to  con¬ 
struct  an  initial  field  for  the  PWE  for  this  problem.  Since  the 
agreement  is  good  for  r„  =  50  m,  the  TDPE  can  handle  this 
problem  using  ph  as  an  initial  condition  at  rn  =  50  m. 


IV.  DISPERSIVE  SEDIMENTS 


Dispersion  and  attenuation  are  introduced  into  the 
TDPE  as  a  perturbation  with  the  definition  K  =  k  +  b, 
where  b4,k  is  a  complex  function  of  to.  To  obtain  a  TDPE 
with  attenuation  and  dispersion,  the  expression  K  =  k  +  b  is 
substituted  into  Eq.  (23)  and  the  Fourier  transform  is  in¬ 
verted.  This  approach  allows  an  arbitrary  complex  disper¬ 
sion  relation. 

Experimental  studies  involving  various  materials  and 
frequencies  have  determined  that  Im (b)  depends  linearly  on 
co  (Refs.  10-12).  In  a  theoretical  study,13  Futterman  showed 
that  R e{b)  is  determined  unambiguously  from  Im(i)  by  the 
principle  of  causality  and  that  the  following  relation  between 
phase  velocity  cp  and  frequency  holds: 


1  1  lap  .  co 

—  - - —  log  — 

Cp  c  ire  co(l 


(57) 


RANGE(kxn) 


I  - 1  ( i  (i  Wiik'-angle  Pf.  results  fur  deep  source.  Solid  curve  ohtamcil  with  l\ 
,tt  (j1  r„  2*0  in  :md  ib,  r,t  50  in.  Diislied  curve  obtained  \v  till  (iimsM.m 
l*K  tu'cr 

Agoi.j1-’  S<-C  Am  .  Vo!  'sin  Hocyi 


where  to,,  is  a  very  low  reference  frequency  and  cp  (co„)  =  c. 
Wuenschel  showed  that  the  agreement  between  prediction 
and  observation  is  excellent  if  both  R e(b)  and  Im(b)  are 
imposed  and  that  the  agreement  is  poor  if  only  Im (b)  is 
imposed."  The  effect  of  sediment  dispersion  on  nearfield 
acoustic  propagation  in  the  ocean  has  been  studied  experi¬ 
mentally,22  and  distortions  attributed  to  sediment  dispersion 
were  observed  in  signals  received  in  the  sediment.  In  this 
article,  we  are  interested  in  the  farfield  effects  of  sediment 
dispersion  on  signals  received  in  the  water  column. 

Adding  the  dispersive  term  from  Eq.  (57)  to  Eq.  (23), 
we  obtain 


dr  \c  c()J  c„ 


liaPto 


rrc„ 


X  log 


u+_2uoctrL_  v 


(58) 


4 co~  -f-  c, |  Z. 

Inverting  the  Fourier  transform  in  Eq.  (58),  we  obtain  the 
wide-angle  TDPE  with  dispersion: 


du 

2c„(d  /dt)L 

lc„  e)  dt  + 

4  (d2/dt2)  —  cf,L 

(59) 

Mu  = 

iJJ  J 

' .  2  ico  .  co  3 

|ru|  -f - log  - 

X  TT  CO„  ) 

Xexp  [ico(i‘ —  t)]dt' dco,  (60) 

where  M  is  the  attenuation/dispersion  operator.  The  equa¬ 
tion 


(61) 


is  solved  numerically  with  an  approach  similar  to  the  ap¬ 
proach  used  to  solve  Eq.  (31 ).  For  a  single  frequency.  Eq. 
(61 )  has  the  exact  solution 


U(r  +  &r,z,co) 


=  U(r,z,co)e\p 


(62) 


where  \  =  cr/?Ar/c„.  The  solution  of  Eq.  ( 6 1 )  is  obtained  by 
inverting  the  Fourier  transform  in  Eq.  (62)  and  substituting 
Eq.  (20)  for  U  to  obtain 


u(r  +  \r,z.t) 


u(r,z,t')e\ p(  -  \'co) 


>  cos 


<o(t'  -  !)  log 


tit '  cho. 
(63) 


The  integral  m  Eq.  (35)  is  approximated  by  dividing  the 
time  integral  into  a  sum  of  integrals  over  small  intervals  over 
which  u  is  assumed  constant.  This  approach  does  not  give  a 
robust  numerical  solution  of  Eq.  (63).  Since  Eq.  (63)  ac¬ 
counts  for  both  dispersion  and  attenuation,  it  must  translate 
and  dampen  a  waveform.  Assuming  u  constant  over  each 
subinterval  of  time  amounts  to  assuming  that  u  is  a  step 
function  in  time.  Thus  one  might  ^Apect  this  approach  to  be 
ineffective  because  the  values  at  the  grid  points  do  not 
change  as  a  step  fuiu.tion  propagates  a  small  distance  A 
robust  sohition  is  obtained  by  approximating  u  by  a  parabola 
over  each  subintei  val 


l.d  1 9HR 
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For  t,„  —  Ar/2  <  t  '</,„+  A/  /2,  we  approximate 
u(r,z,t ')  with 


u(r,z.t')  =  u(r,z,t,„  )  + 


u(r,z,t,„  ,  ,)  -  u(r,z,t,„  ,) 


X  (/'  —  /„,) 

u(r,z,t,„  ,  ,  )  -  2 u(r,z,tm )  +  u(r.z,t,„  ,  ) 

+  2(Af): 

X  (/  ’  —  fm  )2-  (64) 

Substituting  Eq.  (64)  in  Eq.  (63),  we  obtain  the  approxi¬ 
mate  solution 

u(r+  A  r,z.;„) 


V  F(f,„  —  /„  )t/(r.z,/„, ) 


0  too  200  300 


FREQUENCY(Hz) 

HG.  7.  Solid  curve  is  phase  velocity.  Dashed  curve  is  croup  velocity  The 
difference  between  the  curves  is  nearly  constant. 


+  (7(  f„>  ~  6, ) 


u(r.z.t,„  t  ,  )  -  u(r.z,t,„  ,) 


+  H(t,„  -  i„ ) 

u(r.z.l,„  .  |  )  2 u(r.z.t„, )  +  u(r.z.t,„  ,  ) 


F(t)=  —  J  —  exp(  -  paGsin  tit  At  j 

x  cos^ cot  — —  to  log  ]  dco,  ( 66 ) 

Cr(/)=-*-^  —  exp(  -  yco)  At  cos  ^  o)  At  j 

— —  sin  f—  co  At  )  sin  (cot  —  to  log  —  dco. 
co  l  2  /  J  V  tr  coj 

(67) 

H(t)=—(  —  exp(  —  yco)  1—L  sin  (—  co  At  ) 
n  Jo  co  12  \  2  / 

2At  /I  ,  \ 

-I - cos  to  At 

co  V  2  ) 

— sin  o>  At  1 1  cos  (rot  — —  co  log  — 1  dco. 
co~  \2  /  J  V  tr  at,/ 


By  manipulating  the  indices,  Eq.  (65)  becomes 


i/(r+ Ar,z,t„)  =  £  /f,„w(r,2,/m  +  „),  (69) 

m  m  | 

where  the  coefficients  /f,„  are  easily  solved  for  in  terms  of  F, 
<7,  and  //.  Since  the  operator  converges  to  6(  t '  —  t )  as  Ar — 0, 
the  sum  can  be  collapsed  to  a  small  number  of  grid  points 
near  t„  making  the  operator  efficient.  The  integrals  for  F,  G, 
and  H  must  be  evaluated  numerically,  but  this  is  a  relatively 
minor  task. 

To  demonstrate  the  robustness  of  Eq.  (69),  we  allow  it 
to  act  on  a  plane  wave  with  the  initial  condition 

u(r)),z,t )  =  sin(&i|t)  +  3  sin(<u2t),  (70) 

whereat,  =  lOOrrs  'and  at,  =  340rrs‘  '.  The  exact  solution 
is  given  by 


u(r„  +  A  r.z.t) 

=  exp(  —  yco,  )sin  (co,t  —  — i-  aqlog 

V  TT  C0(]  I 

4-  3  exp(  -  paMsin  (co.t  -  —  to.  log  .  (71  ) 

We  let  c  —  1520  m/s.  r„  =  1500  m/s.  ft  =  0.5  dB//., 
to,,  =  2tr  s  Ar  —  2  m.  and  A t  =  1/3  ms.  The  phase  ar.d 
group  velocities  appear  in  Fig.  7.  It  is  easy  to  show  frorp  Eq. 
( 57 )  that  the  difference  between  the  phase  and  group  veloc¬ 
ities  is  nearly  constant  over  a  wide  frequency  band.  The  con¬ 
stants  A„,  appear  in  Table  III  for  jm|<c  20.  In  contrast  to  the 
constants  for  the  attenuation  operator.  A„,  is  negligible  for 
m  >  l  due  to  the  causality  of  the  attenuation/dispersion  op¬ 
erator.  We  take  m,  =  -  20  and  m.  =  1.  The  results  of  Eq. 
(69)  appear  in  Fig.  8  after  100  and  200  range  steps.  The 
solution  in  the  absence  of  dispersion  is  included  for  empha- 


TABI.E  III.  Coefficients  for  the  numerical  solution  of  the  attenuation/dis¬ 
persion  operator  The  small  numbers  in  the  right  column  are  not  used  in  the 
calculations. 


m 

Am 

m 

Am 

-  20 

5.5357370E-05 

2 

2.043 1539E-09 

-  19 

6.395 1957E-05 

3 

3.0856633E-09 

-  18 

7.I207804E-05 

4 

4. 1212398E-09 

-  17 

7.9770936E-05 

5 

5  1469717E-09 

16 

8.9975285E-0S 

6 

6. 1596332E-09 

-  15 

1. 02268  50E-04 

7 

7. 1 556476E-09 

-  14 

1 . 1 726 1 87E-04 

8 

8.1310549E-09 

-  13 

1.3580774E-04 

9 

9.0814698E-09 

-  12 

1.59I2310E-04 

10 

1  0002032E-08 

-  11 

'  .8  899402  E-04 

11 

1. 088737  IE-08 

-  10 

2.2812726E-04 

12 

1 . 1731 533E-08 

-  9 

2.8078392E-04 

13 

1  2527930E-08 

-  8 

3.5398538E-04 

14 

I.3269264E-08 

-  7 

4.5996779E-04 

15 

1.3947457E-08 

-  6 

6  2I60386E-04 

16 

1  4553552E-08 

-  5 

8.8567595E-04 

17 

1.5077628E-08 

-  4 

1.3598768E-03 

18 

1  5508673E-08 

-  3 

2.332 1055E-03 

19 

1. 583447  3F-08 

-  2 

4.5252889E-03 

20 

1  5370508E-08 

-  1 

-  3.9909109E-02 

0 

0.9272352 

1 

9.9426970E-02 
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(b> 


\ 


0  20  40  60 

RELATIVE  TlME(ms) 

FIG.  $  Benchmark  of  allenuaiion/dispersion  numerical  solution.  Solid 
curve  is  the  numerical  solution  Dashed  curve  is  the  exact  solution  Dotted 
curve  is  the  exact  solution  without  dispersion  term  Re(/>).  after  la)  UK) 
range  steps  and  I  h)  200  range  steps 


TABLE  IV.  Data  for  the  dispersion  problem.  CPI'  for  TDPE  calculation 


z„  =  28  m 

z,  —  28  m 

c,,  =  1  SIX)  m/s 

cu  =  1 500  m/s 

p„  —  1  g/cm* 

P,  -o 

= <•„(">•*) 

ph  =  1.5  g/cm' 

ph  =  0.5  dB/4 

d  =  30  m 

>•  -  300  s  1 

r„  —  20  m 

Ar  =  2  m 

Az  =  1  m 

A  /  —  1/3  ms 

<ou  =  2  rr  s  1 

zif  =  100  m 

CPU  =  3.9  h 

sis.  We  observe  thal  the  shape  of  the  dispersed  wave  is  dis¬ 
torted.  The  agreement  between  the  exact  and  numerical  so¬ 
lutions  is  good. 

To  demonstrate  the  effect  of  sediment  dispersion  on 
pulses,  we  consider  the  problem  described  in  Table  IV.  The 
base  sound  speed  is  c(z)  =  ( 1520  +  5 z)  m/s  for  z  -  d  <  10 
m,  c(z)  =  1570  m/s  for  z  —  d  >  10  m  with  <y„  =  2ns  It  is 
evident  from  Fig.  7  that  the  phase  velocity  cp  is  about  40  m/s 
greater  than  c  in  a  wide  frequency  band  centered  about  100 
Hz.  For  comparison,  we  consider  the  case  with  the  frequen¬ 
cy-independent  sediment  sound  speed  c(z)  +  40  m/s.  a 
sound-speed  profile  that  was  assumed  in  an  experimental 
study.’’  The  source  is  Gaussian  with  v  =  300  s“  Contour 
plots  for  the  solution  of  Eqs.  (24)  and  (59)  appear  in  Fig.  9. 
We  see  from  the  waveforms  in  Fig.  10,  which  is  qualitatively 
similar  to  Fig.  3  of  Wuenschefs  paper,  that  the  signal  is  sig¬ 
nificantly  affected  by  the  dispersion. 

It  was  noted  in  Ref.  23  that  the  linear  loss  law  without 
disperison  does  not  provide  good  predictions  for  an  environ¬ 
ment  similar  to  the  one  we  have  just  considered,  and  the 
linear  loss  law  was  challenged.  The  results  in  Ref  23  are 


u: 


75  4 


100  j - 1 - 1 - - - 1 

1890  1985  1980  2025  7070 

RANGE(m) 


Q 

75- 


100  1 - - - 1 

9890  9995  9980  4025  1070 

RANGE(fii) 


A2 


100  1 - T - - - - 

3890  393f>  3980  4025  4070 

RANGE(m) 


FIG.  9.  Snapshots  of  (he  pressure  field  obtained  with  wide-angle  IDPL  wilh  dispersion  near  la)  r  2  kin  and  (hi  r  4  km.  ami  with  wide-angle  1  DPF. 
without  dispersion  near  (c)  r  2  km  and  Id)  r  4  km 
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FIG  10  Waveform  delected  by  receiver  at  z  =  28  m  and  (a)  r  =  2  km  and 
(b)  r  =  4  km  Solid  curve  is  from  wide-angle  TDPE  with  dispersion. 
Dashed  curve  is  from  wide-angle  TDPE  without  dispersion 


based  on  an  analysis  that  neglects  sediment  dispersion. 
Based  on  the  results  of  Wuenschel  and  the  fact  that  sediment 
dispersion  affects  water-borne  signals  significantly  as  we 
have  demonstrated,  it  appears  that  dispersion  might  be  a 
better  explanation  of  the  observations  than  a  new  loss  law. 


V.  CONCLUSIONS 

Being  the  inverse  Fourier  transform  of  the  PE,  the  wide- 
angle  TDPE  is  analytically  equivalent  to  the  wide-angle  PE. 
Like  the  wide-angle  PE,  the  wide-angle  TDPE  is  accurate 
and  efficient  and  easily  handles  range-dependent  propaga¬ 
tion.  Since  the  TDPE  is  an  initial  value  problem  in  range,  it  is 
easy  to  construct  an  accurate  initial  field  for  the  TDPE.  The 
TDPE  handles  sediment  dispersion,  which  can  have  a  signif¬ 
icant  effect  on  water-borne  signals.  The  TDPE  can  be  initia¬ 
lized  with  the  homogenous  half-space  field,  which  is  accu¬ 
rate  and  easy  to  construct.  The  numerical  solution  of  the 
wide-angle  TDPE  involves  the  method  of  alternating  direc¬ 
tions  as  well  as  several  other  .-.tandard  numerical  methods. 
The  refraction  operator  is  solved  with  the  Lax-Wendroff 
scheme.  The  wide-angle  diffraction/density  operator  is 
solved  with  the  finite  element  method,  centered  differences, 
and  Crank-Nieolson  integration.  The  attenuation/disper¬ 
sion  operator  is  solved  with  quadrature  to  evaluate  integrals. 
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APPENDIX:  DEPTH  DISCRETIZATION  WITH 
GALERKIN’S  METHOD 

Galerkin’s  method  with  linear  test  functions  is  applied 
to  discretize  depth  dependence.  We  use  the  approach  used  in 
Ref.  8  with  some  improvements.  The  depth  grid  points  are 
defined  to  be  z,  =  i  A z,  where  A z  is  the  depth  increment.  The 
basis  functions  (z)  vanish  for  |z  —  z,  j  >  A z,  increase  lin¬ 
early  from  0  to  1  over  z,  ,  <  z  <  z, ,  and  decrease  from  1  to  0 
over  z,  <z<z,  +  , .  The  basis  functions  can  be  used  to  ap¬ 
proximate  a  function  by  a  piecewise  linear  function  with  ex¬ 
act  agreement  at  the  grid  points.  To  solve  the  wide-angle  PE, 
we  define  U, \r)  =  U(r,z(),  ©  =*  log(p),  and  0,  =  0(z, )  to 
obtain 


U(r,z )*£  fr, ■(/•)¥,. (z). 

(Al) 

©(zJa^O.'P.fz), 

(A2) 

Galerkin's  method  is  used  to  discretize  depth  dependence  in 
Eq.  (17)  by  requiring  that  the  following  hold  for  all  /: 


d'U  '■ 
drdz2 


p  dz  dr  dz  dr  p  dz  dz  ) 


dz  =  0. 


(A3) 


Equation  (18)  is  obtained  by  substituting  Eqs.  (Al)  and 
(  A2)  into  Eq.  (A3).  The  entries  of  the  matrices  R  and  Sare 
determined  by  the  following  approximations  obtained  with 
Galerkin’s  method: 


a-v  v. .  i  -  2t/,  +  v,  , 

dz  -•  ( Az)J 

^  v.  .  ,  +  417,  +  U,  , 

6 

<70  dU 
dz  dz  --  ■- 

(20,  -  0,  .  t  -  0,  ,  )U, 

2(Az): 

(0,  4  (0,  ,  —©,)«/,  , 

2(Az): 


( A4) 
( A5) 


( A6) 


The  entries  of  the  matrices  A.  B,  and  C  are  also  derived  from 
the  above  difference  formulas. 
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